function [pROI pWSS ti] = polarROI(pMAP, ROI, stat, radialres)

% polarROI takes a polar map of a shear field (pMAP) and the analogous
% cartesian ROI coodinates and remaps the ROI to polar coordinates using
% stats to center the vessel centroid. Improfile is then used to pull 
% the WSS values (pWSS) at the lumen boundaries.

%% calculate ROI in polar coordinates

x = find(ROI(:,2) > 0);
y = find(ROI(:,1) > 0);  % find all nonzero values of x the ROI coords

[t r] = cart2pol(ROI(y,1)-stat.Centroid(1),ROI(x,2)-stat.Centroid(2));

x       = diff(t);
sd      = std(x);
indwrap = find(x<(-3*sd)|x>(3*sd));
x(indwrap,:) = [] ;
if mean(x) < 0         % flip if going the wrong direction
    t = flipud(t);
    r = flipud(r);
end

t     = t * 180/pi;         % transform to degrees and shift 180 degrees
r     = r/radialres + 1;              % shift r +1 b/c plotting in image coords (no zeros)

indstart = find(t>0,1,'first');  % find where t needs to be wrapped... can use indwrap... but will fix later...
indend   = find(t<0,1,'last');   % find if any positive vals at end (ie ~=size)
if indstart > 1
    shift = size(t,1) - indstart + 1;
    t     = circshift(t,shift);
    r     = circshift(r,shift);
else
    shift = size(t,1) - indend;
    t     = circshift(t,shift);
    r     = circshift(r,shift);
end
negtheta = find(t<0);
t(negtheta) = t(negtheta) + 360;

% interpolate values for 0 and 360 degrees
rise   = r(1) - r(end);
run    = (360 + t(1)) - t(end);
slope  = rise/run;
rat2pi = r(end) + slope*(360 - t(end)); % interp value of r@360,0 degrees
t      = [0;t;360];
r      = [rat2pi;r;rat2pi];

%% calculate WSS along pROI

[ti,ri,pWSS] = improfile(pMAP,t,r);
pROI         = [t r];